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Abstract 

The multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) theory within second 
quantization representation of the Fock space, a novel numerically exact methodology to treat 
many-body quantum dynamics for systems containing identical particles, is applied to study the 
effect of vibrational motion on electron transport in a generic model for single-molecule junctions. 
The results demonstrate the importance of electronic-vibrational coupling for the transport char- 
acteristics. For situations where the energy of the bridge state is located close to the Fermi energy, 
the simulations show the time-dependent formation of a polaron state that results in a pronounced 
suppression of the current corresponding to the phenomenon of phonon blockade. We show that 
this phenomenon cannot be explained solely by the polaron shift of the energy but requires meth- 
ods that incorporate the dynamical effect of the vibrations on the transport. The accurate results 
obtained with the ML-MCTDH in this parameter regime are compared to results of nonequilibrium 
Green's function (NEGF) theory. 



I. INTRODUCTION 



Charge transport in single-molecule junctions, i.e. single-molecules that are bound to 
metal or semiconductor electrodes, has been of great interest recently.-"— Employing differ- 
ent experimental techniques, including electromigration or mechanically controllable break 
junctions or scanning tunneling microscopy,-*^"— the conductance properties of nanoscale 
molecular junctions have been investigated. The observed current-voltage characteristics 
typically exhibit a nonlinear behavior with resonance structures at larger bias voltages as- 
sociated with the discrete energy levels of the molecular bridge. The experiments have also 
revealed a wealth of interesting transport phenomena including Coulomb blockade,— the 
Kondo effect,— negative differential resistance,— i^*^ switching and hysteresis.—"— Further- 
more the possibility to obtain transport characteristics that resemble those of a diode^^ or 
a transistor— has been demonstrated. These findings have stimulated great interest in the 
basic mechanisms which govern quantum transport at the molecular scale. 

An interesting aspect that distinguishes single-molecule junctions from mesoscopic de- 
vices is the influence of nuclear motion on electron transport. Because of the small size of 
molecules, the charging of the molecular bridge is often accompanied by significant changes 
of the nuclear geometry that result in strong coupling between electronic and vibrational 
degrees of freedom. This coupling may give rise to substantial current- induced vibrational 
excitation and thus may cause heating and possible breakage of the molecular junction. 
The signature of nuclear motion has been observed in conduction measurements of a vari- 
ety of molecular junctions,— i^ii^ii^ii^>2^"—>^"— e.g., H2 between platinum electrodes," Cgo 
molecules between gold electrodes,— and copper phthalocyanine^^ on aluminum oxide film. 
Vibrational signatures of molecular bridges have also been observed in inelastic electron 
tunneling spectroscopy.—"— New experimental techniques^"— based, e.g., on Raman spec- 
troscopy, allow the characterization of the nonequilibrium state of the vibrational degrees of 
freedom in a molecular junction. 

The experimental progress has stimulated much interest in the theoretical modeling and 
simulation of vibrationally coupled electron transport in molecular junctions. To this end, 
a variety of theoretical approaches have been developed and employed, including scat- 
tering theory,—"— nonequilibrium Green's function approaches,—"— and master equation 
methods.—!^"— Although much physical insight has been obtained by the application of 



these methods, all these approaches involve significant approximations. For example, NEGF 
methods and master equation approaches are usually based on (self-consistent) perturbation 
theory and/or employ factorization schemes. Scattering theory approaches to vibrationally 
coupled electron transport, on the other hand, neglect vibrational nonequilibrium effects 
and are limited to the treatment of a small number of vibrational degrees of freedom. Fur- 
thermore, a systematic improvement of these approaches to yield numerically exact result, 
though formally possible by, e.g., including higher orders in the perturbation expansion, 
is practically very challenging. These shortcomings have motivated us to develop a sys- 
tematic, numerically exact methodology to study quantum dynamics and quantum trans- 
port including many-body effects, in particular, correlated electronic-nuclear dynamics — 
the multilayer multiconfiguration time- dependent Hartree (ML-MCTDH) theory in second 
quantization representation (SQR).''^ Other efforts along the same direction include the nu- 
merical path integral approach,—"— real-time quantum Monte Carlo simulations,—"^ the 
numerical renormalization group approach,— and the time- dependent density matrix renor- 
malization group. For a comparison and an overview of various different methods in the 
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related problem of nonequilibrium transport with electron-electron interaction see Ref. 

In this paper, we report results of accurate quantum simulations employing the ML- 
MCTDH-SQR theory for a generic model of vibrationally coupled electron transport through 
molecular junctions. The paper is organized as follows. Section |II]outlines the physical model 
and the observables of interest. The ML-MCTDH-SQR theory is described in Section lllli 
Section [TVl presents numerical results for vibrationally coupled electron transport in different 
parameter regimes as well as an analysis of the transport mechanism. Moreover, the validity 
of NEGF theory in the regime of phonon blockade is discussed. Finally, Section IVl concludes. 



II. MODEL AND OBSERVABLES OF INTEREST 



To study vibrationally coupled electron transport we consider a simple generic model 
for a single-molecule junction. It comprises one discrete electronic state at the molecular 
junction, two electronic continua describing the left and the right metal leads, respectively, 
and a distribution of harmonic oscillators that models the vibrational modes of the molecular 
bridge. The Hamiltonian reads 

H = Hel + H^uc + Hel-nuc, (2.1a) 



where Hei, -f^nuc, and Hei-nuc describe the electronic degrees of freedom, the nuclear vibra- 
tions, and their coupling terms, respectively 
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-?^oi-nuc = d'^dY^CjQj. (2. Id) 
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Thereby, d'^/d, c^^/ck^, cl^/ckj^ are the fermionic creation/annihilation operators for the 
electronic states on the molecular bridge, the left and the right leads, respectively. The 
corresponding electronic energies E^^, E^^ and the molecule- lead coupling strengths V^fc^, 
Vdfc^, are defined through the energy- dependent level width functions 

Tl{E) = 271 Y,\ydk,?KE-Ek^), r^(E) = 27r J]|\/rffcj2(5(E-E,^). (2.2) 

In principle, the parameters of the model can be obtained for a specific molecular junction 
employing first-principles electronic structure calculations.^® In this paper, which focuses on 
the transport methodology, however, we will use a generic parametrization. Employing a 
tight-binding model, the function T{E) is given as 

. ( ^ ^l^'l . (2,3a) 
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Tl{E) = T{E - /ii), r^(^) = T{E - fiR), (2.3b) 

where /3e and are nearest-neighbor couplings between two lead sites and between the lead 
and the bridge state, respectively. I.e., the width functions for the left and the right leads 
are obtained by shifting T{E) relative to the chemical potentials of the corresponding leads. 
We consider a simple case of two identical leads, in which the chemical potentials are given 
by 

liL/R = Ef± V/2, (2.4) 

where V is the bias voltage and Ef the Fermi energy of the leads. Since only the difference 
Ed — Ef is physically relevant, we set E'j = in this paper. 



The frequencies Uj and electronic-nuclear coupling constants Cj of the vibrational modes 
of the molecular junctions are modeled by a spectral density function^"^ 



9E77^(^-^^)- (2.5) 



2 UJ. 



In this paper, the spectral density is chosen in Ohmic form with an exponential cutoff 

Mu) = \auje--l-% (2.6) 

where a is the dimensionless Kondo parameter. 

Both the electronic and the vibrational continua can be discretized using an appropri- 
ate scheme.— Within a given time scale the numbers of electronic states and bath modes 
are systematically increased to reach converged results for the quantum dynamics in the 
condensed phase. In this paper, we employ 64-128 states for each electronic lead, implying 
32-64 electrons per lead, and a bath with 100-400 modes. 

The observable of interest in transport through molecular junctions is the current for a 
given bias voltage, given by (in this paper we use atomic units where ^ = e = 1) 

Hit) = = -^tr {pe^^H[H, iV,]e-^*} , (2.7a) 

Init) = = ^tr {pe^^H[H, NnK^^'} • (2.7b) 

Here, Nijjiit) denotes the time-dependent charge in each lead, defined as 

Ndt) = ^^tT[pe'^'N^e-'^% C = L,R. (2.8) 

In the expression above = c^^c^^ is the occupation number operator for the electrons 
in each lead {( = L, R) and p is the initial density matrix representing a grand-canonical 
ensemble for each lead and a certain preparation for the bridge state 



'° exp [-PiHo - PlNl - PrNr)\^ , (2.9a) 



P = Pd 

Ho = Yl ^^A^^^ + E E^rctcu. + ^nuc- (2.9b) 

Here /3|] is the initial reduced density matrix for the bridge state, which is usually chosen as 
a pure state representing an occupied or an empty bridge state, and H'^^^ defines the initial 
bath equilibrium distribution, e.g., ^^nuc given above. The dependence of the steady-state 



current on the initial density matrix is a complex issue and is partly addressed in the results 
section of this paper. Finally, for Hamiltonian (12. ip the explicit expression for the current 
operator is given as 

= z[H, N^] = iJ2 Vdk,{d+Ck^ - 4d), ( = L,R. (2.10) 

The transient behavior of the thus defined currents /i?,(t) and Iiit) is usually different. 
However, the long-time limits of Init) and Iiit), which define the stationary current, are 
the same. It is found that the average current 

m = \[In{t) + lL{t)], (2.11) 

provides better numerical convergence properties by minimizing the transient characteristic, 
and thus will be used in this paper. 

Within the model, the left and right leads represent the electronic continuum, each con- 
taining an infinite number of states. As mentioned above, in our simulations this continuous 
distribution is represented by a finite number of electronic states. The number of states 
required to properly describe the continuum limit depends on the time t. The situation is 
thus similar to that of a quantum reactive scattering calculation in the presence of a scat- 
tering continuum, where, with a finite number of basis functions, an appropriate absorbing 
boundary condition is added to mimic the correct outgoing Green's function.—"— Employing 
the same strategy for the present problem, the regularized electric current is given by 

r^s= hm r dt'^e-'^K (2.12) 

The regularization parameter r] is similar (though not identical) to the formal convergence 
parameter in the definition of the Green's function in terms of the time evolution operator 

POO 

G{E+) = Mm i-i) dte'^^+''^-"^K (2.13) 

In numerical calculations, rj is chosen in a similar way as the absorbing potential used in 
quantum scattering calculations.—"— In particular, the parameter rj has to be large enough to 
accelerate the convergence but still sufficiently small in order not to affect the correct result. 
While in the reactive scattering calculation rj is often chosen to be coordinate dependent, in 
our simulation rj is chosen to be time dependent 

,0 {t<r) 
r^(t) = { ^ ) (2.14) 



Here rjo is a damping constant, r is a cutoff time beyond which a steady state charge ffow 
is approximately reached. As the number of electronic states increases, one may choose a 
weaker damping strength rjo and/or longer cutoff time r. The former approaches zero and 
the latter approaches infinity for an infinite number of states. In practice, for the systems 
considered in this work, convergence can be reached with a reasonable number of electronic 
states in the range of 64-128, with a typical r = 40-60 fs (a smaller r for less number of 
states) and l/?7o = 2-5 fs. 

To analyze the transport mechanisms, it is also expedient to consider the population of 
the electronic state localized on the the molecular bridge, which is given by 



III. THE MULTILAYER MULTICONFIGURATION TIME-DEPENDENT 
HARTREE THEORY IN SECOND QUANTIZATION REPRESENTATION 

The accurate treatment of the time-dependent transport problem outlined above requires 
a method that is able to describe the quantum dynamics of a system with many electronic 
and nuclear degrees of freedom. To this end, we employ the recently proposed Multilayer 
Multiconfiguration Time-Dependent Hartree Theory in Second Quantization Representation 
(ML-MCTDH-SQR).-2£ This method extends the ML-MCTDH approach to the treatment 
of indistinguishable particles. 

A. General formulation of the ML-MCTDH theory 

The ML-MCTDH theory— is a rigorous variational method to propagate wave packets 
in complex systems with many degrees of freedom. In this approach the wave function is 
represented by a recursive, layered expansion. 




(2.15) 
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(3.1a) 



Q(k) 



(3.1b) 



M(K,g) 

where Aj-^^j^ j^(t), B^^ff^^^^^^^^it), Cafa2-..aM{n.q){'^) so on are the expansion coefficients for 
the ffist, second, third, layers, respectively; \ip'^^{t)), \v\^''^\t)), {t)) , are the 

"single particle" functions (SPFs) for the ffist, second, third, layers. In Eq. flS.lal) . p 
denotes the number of single particle (SP) groups/subspaces for the ffist layer. Similarly, 
Q{k,) in Eq. fl3.1bp is the number of SP groups for the second layer that belongs to the 
Kth SP group in the ffist layer, i.e., there are a total of X]k=i Q(^) second layer SP groups. 
Continuing along the multilayer hierarchy, M{k,, q) in Eq. (13. Id) is the number of SP groups 
for the third layer that belongs to the gth SP group of the second layer and the Kth SP group 
of the ffist layer, resulting in a total of Y.l=i Hqif MiK, q) third layer SP groups. Naturally, 
the size of the system that the ML-MCTDH theory can treat increases with the number of 
layers in the expansion. In principle, such a recursive expansion can be carried out to an 
arbitrary number of layers. The multilayer hierarchy is terminated at a particular level by 
expanding the SPFs in the deepest layer in terms of time-independent configurations, each 
of which may contain several Cartesian degrees of freedom. 

The variational parameters within the ML-MCTDH theoretical framework are dynami- 
cally optimized through the use of Dirac-Frenkel variational principle^I 

{8^>m^-H\^>{t))=^, (3.2) 

which results in a set of coupled, nonlinear differential equations 

Ll coefficients 

H{t)\m{t)), (3.3a) 

#'^(t))L2 coefficients = [1 - P^^\t)][p^^\t)]-\H)^^\t)\^^'^\t)) , (3.3b) 
^|^!(^'^^^))L3 coefficients = [1 " ^S"'^ (^)] [^^'^ W] (^)ir/Hi) 1^^"''^^^) ) , (3.3c) 
^lf ''•'\t))L4coeffi.ents = [1 - (t)] [^fe^'"^ (t)] (H)^"'^ (t) l^^'^'^'^^t) ) , (3.3d) 



For referencing purpose we label the SP subspaces from top to bottom layers as level one 
(Ll), level 2 (L2), and so on. In the notation used in Eq. (13. 3 p the time derivatives on the left 



hand side (represented by a dot) are only performed with respect to the expansion coefficients 
of a particular layer (denoted by the respective subscript). For example, the time derivative 
in Eq. fl3.3ap acts only on the LI expansion coefficient Aj^^j,^ j^it); the time derivative in 
Eq. fl3.3b|) is on the L2 expansion coefficient -Bf^f^" «q( -)^^)' "'"^ "^^^ convention, for a 

A^-layer version of the ML-MCTDH theory there are + 1 levels of expansion coefficients. 
In this sense the conventional wave packet propagation method is a "zero-layer" MCTDH 
approach. 

In practical implementations various intermediate quantities are defined within the sub- 
spaces of each layer of the wave function.— For example, the top-layer Hamiltonian matrix 
[Hit)]jL = ($j(t)|i^|$L(t)), where |$j(t)) = HLi I'/'j^W) is a configuration in the LI 
subspace; the reduced density matrices p^''-'(t), gi^2'^\t), Qi^s''''^\t), Qi^^^'"''"'\t) for the 
ffist, second, third, and Nth layers, respectively; and the corresponding mean- field opera- 
tors (i^)(")(t), {n)L2''\t), Wtt^\t), (^)LN'^""^(i)- These operators can be recursively 
evaluated by means of the single hole functions \^!^m{t)), |5'L2'.^m s(^))' IS'Ls'^mll^))' 
ffist, second, third, and further layers.— i^i^^ 

The introduction of this recursive, dynamically optimized layering scheme in the ML- 
MCTDH wavefunction provides more flexibility in the variational functional, which results 
in tremendous gain in our ability to study large quantum many-body systems. During the 
past few years, significant progress has been made in further development of the theory 
to simulate quantum dynamics and nonlinear spectroscopy of ultrafast electron transfer re- 
actions in condensed phases.— i^^i"— The theory has also been generalized to study heat 
transfer through molecular junctionsi^E and to calculate rate constants for model proton 
transfer reactions in molecules in solution.—"^ Recent work of Manthe has introduced an 
even more adaptive formulation based on a layered correlation discrete variable representa- 
tion (CDVR).^'ii^ 

B. Treating identical particles using the second quantization representation of 
Fock space 

Despite its previous success, the original ML-MCTDH theory was not directly applicable 
to studying systems of identical quantum particles. This is because an ordinary Hartree 
product in the ffist quantized picture is only suitable to describe a configuration for a 



system of distinguishable particles. To handle systems of identical particles, one strategy 
is to employ a properly symmetrized wave function, i.e., permanents in a bosonic case 
or Slater determinants in a fermionic case. This led to the MCTDHF approach^"— for 
treating identical fermions and the MCTDHB approacb^i^ for treating identical bosons as 
well as combinations thereof.-^^ However, this wave function-based symmetrization is only 
applicable to the single layer MCTDH theory but is incompatible with the ML-MCTDH 
theory with more layers — there is no obvious analog of a multilayer Hartree configuration 
if permanents/determinants are used to represent the wave function. As a result, the ability 
to treat much larger quantum systems numerically exactly was severely limited. 

To overcome this limitation we proposed a novel approach^ that follows a fundamentally 
different route to tackle quantum dynamics of indistinguishable particles — an operator- 
based method that employs the second quantization formalism of many-particle quantum 
theory. Thereby the variation is carried out entirely in the abstract Fock space using the oc- 
cupation number representation. This differs from many previous methods where the second 
quantization formalism is only used as a convenient tool to derive intermediate expressions 
for the first quantized form. In the new approach the burden of handling symmetries of 
identical particles in a numerical variational calculation is shifted completely from wave 
functions to the algebraic properties of operators. 

The procedure can be illustrated by considering a system of identical fermions. In the first 
quantized representation a Slater determinant \xpiXp2---XPn) describes an anti-symmetric 
A^-particle state by choosing spin orbitals out of the M orthonormal spin-adapted basis 
functions, {xi(x), X2(x), Xm(x)}. All possible Slater determinants in this form constitute 
the fermionic subspace of the A^-particle Hilbert space, denoted by Tii^M, N). The Fock space 
J-'(M) is formed by considering an arbitrary number of particles 

jr(M) = n{M, 0) © n{M, 1) © n{M, 2) © ... © n{M, M). (3.4) 

In second quantization a convenient basis to represent the Fock space is the occupation 
number basis 

|n) = |ni,n2, ...,nM), (3.5) 

where np can be either 1 if the one-particle state xp is occupied [i.e., present in the original 
Slater determinant] or if it is unoccupied. The number of particles in state |n) is given by 
A^ = ^Py ^-iid thus T-L{M, N) contains all occupation-number vectors with A^ particles. 



Each occupation-number state is defined by acting a series of creation operators ap on tlie 
vacuum state, |n) = Y[p=ii'^p)^''\^^^) ■ 

In contrast to a Slater determinant, tlie occupation- number state |n) can be formally 
written in the form of a Hartree product,— 

M 

|n) = Y\_ I'^p) = \ni)\n2)...\nM)- (3.6) 

p=i 

This formal factorization suggests a different decomposition of the Fock space 

T{M) = ® Ml) ® ... ® Ul) ® ... ® /a/ (1), (3.7a) 

where /k(1) represents a single spin-orbital subspace with two possibilities/states: occupied 
or unoccupied. Bigger subspaces can be formed by grouping a few spin orbitals together 

L 

J-(M) = /i(mi)®/2(m2)®...®/.(m,)®...®/L(mi), M = J]m«, (3.7b) 

K=l 

where fn^fn,^) represents the subspace with spin orbitals and thus 2^'^ states. 

The decomposition schemes in Eqs. f l3.7al) and (I3.7b[) are conceptually different from that 
in (13. 4p . They no longer require dealing with a full-length occupation-number vector in one 
step and treating it as a whole, unbreakable object like the original Slater determinant, 
but rather focus on each subspace /k(^k) containing 2™"= linearly independent sub- vectors 
{0^^''}, In = 1, 2™". The crucial point is that in a variational calculation each sub- vector 
in the nth subspace is not restricted to a particular fixed basis vector {4>\'^J}, but may also be 
an appropriate superposition of all the sub-vectors in this subspace. One may then express 
the overall wave function in the same multilayer form as in Eq. (13. ip . Take, for example, 
three explicit layers in (13. ip . the deepest layer (3rd layer here) is expanded in the full basis 
sub-vectors in the Fock subspace as, 

\Crit)) = E - E ^n;nI'X(«.,.,(^) \ni)\n2)-\n^i.,,,,)). (3.8) 

ni=0n2=0 rj„(«_,,^)=0 

|\l/(t)) in Eq. (13.10 is then built "bottom-up" from the optimal, time-dependent SPFs for 
the subspaces. 

After introducing the occupation-number representation of the Fock space, the Hamil- 
tonian can be expressed in terms of fermionic creation/annihilation operators via the stan- 
dard procedure in the second quantization formalism.—*^ The overall method is thus a 



ML-MCTDH theory in second quantization representation (SQR). The major difference be- 
tween the ML-MCTDH-SQR theory for identical fermions and the previous ML-MCTDH 
theory for distinguishable particles is the way how operators act. In the second quantized 
form the fermionic creation/annihilation operators fulfill the anti-commutation relations 



{ap, Qq} = apOQ + ttgOp = 5pQ, {op, Oq} = {ap, aq} = 0. 



(3.9) 



The symmetry of identical particles is thus realized by enforcing such algebraic properties 
of the operators. 

The practical procedure can be illustrated by considering a single layer theory in the form 
of Eq. (I3.7bl) . where each SP group k corresponds to a Fock subspace in (I3.7bl) 



1 1 



K = l 



IV^S^W) = E = E E - E ^nX-nr^S^) \n,)\n2)...\n^J. (3.10b) 

/k=1 ni=0n2=0 ?imK=0 

Without loss of generality let us consider acting a creation operator (a^T'')"'" on the SPFs. In 
practical implementation this operation is equivalent to 



(3.11) 



where (/i = 1, 2, /t — 1) is the permutation sign operator that accounts for permuting 
(a[r^)^ from the first subspace all the way through to the nth subspace, and (ai^^)^ is the 
reduced creation operator that only takes care of the fermionic anti-commutation relation 
in the nth subspace. The operator-based anti-commutation constraint (13. 9p results in the 
following operations 



1 1 



1 



(4 



ni=0n2=0 nmK=0 
1 1 1 



V-1 



n(-i 



.'2=1 



s>5-:'w> = E E - E 



.9=1 



K'X.n^St) \n^)\n,)...\l,)...\n^J, 

(3.12a) 

i?C;^2...n™„W 1^1)1^2). (3.12b) 



I.e., (Si/ )"•" not only creates a particle in the z/th spin orbital if it is vacant, but also affects the 
sign of each term in this SPF according to where u is located and what the occupations are 
prior to it. Furthermore, the permutation sign operators S^, = 1, 2, k — 1, incorporate 



the sign changes of the remaining spin orbitals in all the SPFs whose subspaces are prior to 
that of (al'^V. 

The implementation of Eq. (13. lip is sophisticated but nevertheless a routine practice in 
the MCTDH or ML-MCTDH theory — a product of operators. Thereby, the action of each 
Hamiltonian term (product of creation/annihilation operators) can be split into a series of 
operations on individual Fock subspaces. 

The generalization from the single layer to the multilayer case is tedious but straight- 
forward. The Fock space is decomposed in a recursive, layered fashion — the spin orbitals 
here in the ML-MCTDH-SQR theory are treated in the same way as the degrees of free- 
dom in the original ML-MCTDH theory, except that the orderings of all the SP groups in 
all layers need to be recorded and maintained in later manipulations. The wave function 
can then be recursively expanded via Eq. (13. ip . More importantly, the equations of mo- 
tion have the same form as in the original ML-MCTDH theory. The only difference is that 
each creation/annihilation operator of the Hamiltonian is effectively a product of operators: 
a reduced creation/annihilation operator that only acts on the bottom-layer SPFs for the 
Fock subspace it belongs to, and a series of permutation sign operators that accounts for 
the fermionic anti-commutation relations of all the spin orbitals prior to it. 

In the second quantized form, the wave function is represented in the abstract Fock 
space employing the occupation number basis. As a result, it can be expanded in the same 
multilayer form as that for systems of distinguishable particles. It is thus possible to extend 
the numerically exact treatment to much larger systems. The symmetry of the wave function 
in the first quantized form is shifted to the operator algebra in the second quantized form. 
The key point is that, for both phenomenological models and more fundamental theories, 
there are only a limited number of combination of fundamental operators. For example, in 
electronic structure theory only one- and two-electron operators are present. This means 
that one never needs to handle all, redundant possibilities of operator combinations as 
offered by the determinant form in the first quantized framework. It is exactly this property 
that provides the flexibility of representing the wave functions in multilayer form and treat 
them accurately and efficiently within the ML-MCTDH-SQR theory. It is also noted that 
the ML-MCTDH-SQR approach outlined above for fermions has also be formulated for 
bosons or combinations of fermions, bosons and distinguishable particles.— Here, we apply 
it to vibrationally coupled electron transport, which involves a combination of vibrational 



degrees of freedom and indistinguishable electrons. 

IV. RESULTS AND DISCUSSION 

In this section we present applications of the ML-MCTDH-SQR methodology to vibra- 
tionally coupled electron transport employing the model described in Sec. HIl We first con- 
sider a model with the following set of electronic parameters: The energy of the discrete 
state Ed is located 0.5 eV above the Fermi energy of the leads {Ef = 0). The tight-binding 
parameters for the function T{E) are = 0.2 eV, /3e = 1 eV, corresponding to a moderate 
molecule-lead coupling an a bandwidth of 4 eV. 

Figure [U shows the time-dependent current for low bias voltage, V = 0.2 V, and a range 
of different temperatures, - 300 K. Panel (a) depicts the purely electronic current obtained 
without coupling to the vibrational degrees of freedom {a = 0). In this case significant 
electronic coherence is observed for the current I{t) at short time, which decreases for longer 
time. A plateau of I{t) is reached in relative short time (~ 30 fs), which demonstrates the 
feasibility of using a time-dependent approach to obtain the stationary current. In contrast 
to the transient characteristics, for the purely electronic problem considered in Fig. [TJ^a) the 
stationary current can also be obtained exactly from the Green's function or the scattering 
theory approach, employing e.g. Landauer theory.- '^^^i-*^^^ The thus obtained stationary value 
of the current agrees with the simulation result. The results in Fig. [T](a) also illustrate the 
fairly weak temperature dependence of I{t) for this set of electronic parameters. 

Figure [H^b) shows results for the same electronic parameters as in Fig. [H^a), however 
including the coupling to the vibrational bath. The characteristic frequency of the bath 
has been chosen as Uc = 500 cm~^ and the overall electronic-nuclear coupling strength is 
determined by the reorganization energy, A = 2auJc = 2000 cm~^. These parameters repre- 
sent typical values for polyatomic molecules. It is noted that in contrast to many previous 
treatments of vibrationally coupled electron transport in molecular junctions, the present 
method allows a nonperturbative, in principle numerically exact treatment of this nonequi- 
librium problem. The results show that the inclusion of electronic-vibrational coupling has 
a significant effect on the transport characteristics. In particular, it causes a quenching 
of the electronic coherence. As a result, the time scale on which the current I{t) reaches 
its stationary value is shorter. Furthermore, the temperature dependence of the current is 



more pronounced than in the corresponding purely electronic case. It is also noted that in 
this particular physical regime the value of the current is larger than for purely electronic 
transport {vide infra). 

Figure [2] shows the time- dependent current for different electronic-vibrational coupling 
strengths. While smaller electronic-vibrational coupling (A < 1000 cm"^) causes mostly 
decoherence in the transient regime, larger coupling is seen to influence also the stationary 
value of the current significantly. This can be understood qualitatively from the fact that the 
coupling to the vibrations effectively lowers the energy level of the bridge (polaron shift). 
For the given voltage, the bare energy of the electronic bridge state is still outside the 
conductance window, which is defined by the chemical potentials of the two electrodes. The 
coupling to the vibrations brings the level closer to the chemical potential of the electrodes. 
As a result, the current is enhanced. The value of this polaron shift of the energy is given 
by the reorganization energy A. For example, while for a value of A = 1000 cm~^ the 
predominant transport mechanism is nonresonant tunneling, for a value of of A = 4000 
cm~^ the polaron-shifted bridge state is already inside the conductance window between the 
chemical potential of the electrodes and thus the transport mechanism is resonant tunneling. 

We next consider a model with the same parameters except that the energy of the bridge 
state is below the Fermi energy, Ed — Ef = —0.5 eV. As shown in Figure |3l in this case 
an increase in the electronic-vibrational coupling strength not only quenches the electronic 
coherence but also reduces the current monotonically. This is due to the fact that if the 
bridge state is located below the Fermi levels of the leads (Figure [3]) the coupling to the 
vibrational bath will shift its effective energy further away from the resonant transport 
regime. Note that, as implied by the model Hamiltonian, Eq. (12. ip . the calculations depicted 
in Figure [3] use the same electronic reference state as for Figure Ej i.e. the polaron shift is 
to lower energies for both calculations. Figure H] displays the current- voltage characteristics 
for the two models with the same electronic parameters as in Figs. [3] and HI respectively, 
coupled to a bath with a characteristic frequency of Uc = 500 cm~^ and a reorganization 
energy of A = 2000 cm^^. The simulation results are compared with the results for a purely 
electronic system obtained using the Landauer formula. The results show a pronounced 
influence of the vibrational coupling. In this particular parameter regime, for Ed — Ef = 0.5 
eV the vibrationally coupled transport current is higher than that for the purely electronic 
model, whereas for Ed — Ef = —0.5 eV the situation is opposite. Although this can be 



qualitatively explained by the polaron shift of the energy of the bridge state as discussed 
above, the actual quantitative prediction is more complex and requires a non-perturbative, 
accurate approach {vide infra). 

It is worthwhile to point out that the initial condition used in the expression for the 
current, Eq. ( 12.9^ . is not unique. For example, one may choose an initially unoccupied 
bridge state and an unshifted bath of oscillators, i.e. ifnuc as given in Eq. (12.11) . On the 
other hand, one may also start with a fully occupied bridge state and a bath of oscillators 
in equilibrium with the occupied bridge state 



Other initial states may also be prepared. Thus the question arises, whether the stationary 
current depends on the initial state that is used in the time-dependent simulation. 

For the model parameters considered in this paper, our calculations show that the sta- 
tionary state is independent on the initial condition. As an example. Figure |5]^a) shows 
that the two initial states discussed above indeed give the same stationary current for the 
electronic parameters = 0.2 eV, /3e = I eV, — Ef = 0.5 eV, = 0.1 V, and the 
vibrational parameters A = 2000 cm~^ and Uc = 500 cm~^, despite the fact that their initial 
transient characteristics are quite different. As illustrated in Figure Mjo), this is due to 
the fact that although the initial bridge state populations Pd(t), defined by Eq. (I2.15p . are 
quite different, they attain the same stationary value within a relatively short time scale. 
For other parameters, however, the stationary state may depend on the initial condition. 
The investigation of the corresponding phenomenon of bistability^"— will be the subject 
of future work.— 

We note that different sets of initial conditions also affect the time scale at which the 
current I{t) reaches its stationary value, as is evident from Figure [5l In our simulations we 
typically choose initial conditions that are close to the final steady state, e.g., an unoccupied 
initial bridge state with an unshifted bath of oscillators for the transport calculations of 
Fig. m^a) and an occupied bridge state with a bath of oscillators in equilibrium with the 
occupied bridge state for calculations of Fig. Mjo). 

We finally consider a model where the energy of the bridge state is located at the Fermi 
energy of the leads. This parameter regime is particularly interesting, because already for 
small bias voltage the transport mechanism corresponds to resonant tunneling and involves 




(4.1) 



mixed electron/hole transport. For a purely electronic model, the Landauer formula predicts 
the maximum current when Ed = Ej. Including the couplings to the vibrational modes, 
however, may have a significant impact on the electric current. This is illustrated in Figure [6] 
for different electronic-vibrational coupling strengths A. It is seen that for short time the 
current I{t) obtained for finite A follows the current for the purely electronic model (A = 0). 
However, after a short transient time the coupling to the vibrations becomes effective and 
results in a suppression of the current. In particular for larger vibrational coupling, A = 
2000 — 4000 cm~^, the effect is very pronounced and the stationary current is essentially 
blocked over a significant range of bias voltages, as is demonstrated by the current-voltage 
characteristics in Figure [71 

The underlying mechanism can be qualitatively rationalized by considering the energy 
level of the bridge state. For any finite bias voltage, the bare energy of the bridge state 
{Ed — Ef = 0) is located between the chemical potential of the leads and thus, within 
a purely electronic model, current can fiow. The coupling to the vibrations results in a 
polaron shift of the energy of the bridge state. For electronic-vibrational coupling strengths 
A > \V\/2 the polaron-shifted energy of the bridge state is below the chemical potentials of 
both leads and thus current is blocked. This effect, referred to as phonon blockade of the 
current, has been observed e.g. in quantum dots.— 

Although the interpretation of the phonon blockade in terms of the energetics of the 
bridge state is appeahng, it should be emphasized that the mechanism of phonon blockade 
involves the formation of a many-body polaron-type state that is significantly more complex 
than this purely electronic picture and cannot be fully described by just considering the 
static shift of the energy of the bridge state. The bare energy and the polaron-shifted 
energy of the bridge state are only two special points on the multidimensional potential 
energy surface of the charged state given by Vrf(Q) = Ed + ^ J2j '^'jQ'j + J2j ^jQj- For values 
A > \V\/2 the potential energy surface of the discrete state crosses the chemical potential 
of the leads as a function of the nuclear coordinates Qj. In this parameter regime, an 
accurate description of the vibrational dynamics and its coupling to the electronic degrees 
of freedom is required to obtain a quantitative description of the many-body polaron state 
and its transport characteristics. This is demonstrated in Figure [HI which compares the 
electric current obtained with a full vibrationally-coupled many-body ML-MCTDH-SQR 
calculation to that obtained with a purely electronic model for a polaron-shifted energy 



of the bridge state, i.e. ^ — A. The comparison shows that the effect cannot be 
described properly with a purely electronic model but requires methods that incorporate 
the dynamical effect of the vibrations on the transport. 

As discussed in the introduction, a variety of approximate methods have been developed 
and employed to describe vibrationally coupled electron transport in molecular junctions, 
including scattering theory,—"— nonequilibrium Green's function (NEGF) approaches,—"— 
and master equation methods.—"^"— However, only very few of them are applicable to the 
present model. This is because the model involves vibronic coupling to a relatively large 
number of vibrational modes in nonequilibrium. Master equation methods are limited to a 
small number of vibrational degrees of freedom that are treated in full nonequilibrium, while 
many scattering theory approaches neglect vibrational nonequilibrium effects. NEGF theory 
is, in principle, applicable to describe coupling to a larger number of harmonic vibrational 
modes in nonequilibrium. However, if implemented within the self-consistent Born approxi- 
mation it is limited to small electron- vibrational coupling and has so far mostly been applied 
in the nonresonant tunneling regime.— To treat vibrationally coupled electron transport in 
the resonant regime, another NEGF method has been proposed by Galperin et al.^^ and 
extended by Hartle et al.— i^*^ Being based on the polaron transformation, this NEGF 
method is in principle able to treat moderate vibronic coupling strengths. We have applied 
this method to the present model employing 10 vibrational modes to model the vibrational 
distribution described by the spectral density, Eq. (12. 6 p .— Fig. |9] shows a comparison of re- 
sults of the NEGF method for the stationary current with those obtained from ML-MCTDH 
for different vibronic coupling strength. Overall, the results indicate that the NEGF method 
is capable of describing the suppression of the current due to phonon blockade. For small 
vibronic coupling the NEGF results are in almost quantitative agreement with the numeri- 
cally exact results. However, for larger vibronic coupling NEGF theory underestimates the 
effect of phonon blockade. This is presumably due to the fact that the present model does 
not exhibit a strict time-scale separation between the electronic and nuclear degrees of free- 
dom, which is a prerequisite for the NEGF method. These results demonstrate that with 
approximate methods the simulation of transport properties of the present model, which 
is in the nonperturbative regime, is challenging. A more general validation of the NEGF 
method in a broader regime requires extensive studies by both the ML-MCTDH-SQR and 
the NEGF methods and will be the subject of future work. 



While the effect of phonon blockade of the stationary current is to be expected for ener- 
getic reasons, the treatment with the ML-MCTDH method also allows a detailed the study 
of the time-dependent formation of the underlying many-body polaron state. For example, 
Figure [10] shows that the transient dynamics depends significantly on the characteristic fre- 
quency of the vibrational bath, Uc, whereas the value of the stationary current is relatively 
insensitive to Uc- This is due to the fact that the frequency Uc determines the timescale on 
which the system moves on the potential energy surface from the initially prepared state, 
corresponding to the bare energy of the bridge states, to the relaxed state below the chemical 
potential of the leads. 

It is worthwhile to emphasize that the mechanism of the phonon blockade analyzed here 
is different from that of the previously discussed Franck-Condon blockade,—"— which also 
leads to a suppression of the current due to strong electronic-vibrational coupling. While 
the former can be removed by a gate potential that shifts the energy of the bridge state into 
the conductance window, i.e. between the chemical potential of the electrodes, the latter is 
rather insensitive to a gate potential. 

We finally discuss some technical details of the numerical calculations employing the ML- 
MCTDH-SQR method. For the parameter regimes investigated in this paper, the stationary 
current is usually reached at approximately 20-60 fs. To ensure convergence most calcula- 
tions were carried out to 100 fs. Within this time scale of simulation, 64-128 discrete states 
are used to represent each lead's electronic continuum, resulting in a total number of 64-128 
electrons in the ML-MCTDH-SQR numerical treatment. The nuclear bath is represented by 
100-400 modes. The converged number of basis functions for these vibrational modes ranges 
from a few to a few hundred. The calculation was performed with a four-layer ML-MCTDH- 
SQR theory, with one top-layer SP group for the vibrational modes and three top-layer SP 
groups for the electronic part. These SP groups are then recursively expanded via a binary 
tree (i.e., two lower-layer SP groups in each preceding upper- layer SP). The final converged 
results (to within 10% relative error) require 40-60 SPFs for each top-layer electronic SP 
group, 20-40 SPFs for all lower-layer electronic SPs, and 10 SPFs for all nuclear SP groups. 
This results in a total of ~ 10^ equations to solve. Each simulation took between 20 hours 
and several days of CPU time on a typical personal computer. Calculations for a finite 
temperature bath requires ensemble average over a few hundred initial wave functions, and 
were performed on a Cray XT4 parallel computer. 



The convergence of the ML-MCTDH-SQR simulation is illustrated in Figure [TTl where 
the physical parameters are the same as those in Fig. [1] except for a different voltage of O.IV. 
The calculations are performed with a four-layer ML-MCTDH-SQR scheme. For simplicity 
the convergence is shown for the following three different categories. In Figure [TTT a) there 
are 64 electronic states for each lead and 50 modes for the nuclear bath. The number of 
the single particle functions (SPFs) is 10 for each nuclear SP group (of each layer), and 
the number of the SPFs for each electronic SP group (of each layer) is set the same and 
varied. In Figure [TTT b) there are 64 electronic states for each lead, with 40 SPFs for each 
electronic SP group. The number of bath modes is varied, with 10 SPFs for each SP group 
of each layer. In Figure [TTT c) there are 50 bath modes, with 10 SPFs for each SP group. 
The number of electronic states for each lead is varied, with 40 SPFs for each SP group (of 
each layer). 

First, we consider the convergence with respect to the number of SPFs. Fig. ITlT a) shows 
that the electronic (fermionic) degrees of freedom require a relatively large number of SPFs to 
achieve convergence. The steady state current obtained with 20 SPFs differs approximately 
20% from the converged value, although the short time transient dynamics agrees well with 
other results obtained with a larger number of SPFs. In this case convergence is reached when 
the number of SPFs exceeds 30 for the electronic degrees of freedom. For the nuclear degrees 
of freedom, tests have shown that results obtained with 6 or 8 SPFs are nearly identical to 
that obtained with 10 SPFs. Thus, although we used 10 SPFs in all the calculations, we 
believe a smaller number could be equally satisfactory. 

Next, we check convergence with respect to the number of bath modes. Figure [TlT b) 
shows that only small differences are found when the number of modes changes from 10 to 
200. If one is only interested in the steady state current and not the finer details of the 
transient dynamics of I{t), then a bath of 10 modes is sufficient as is used in our NEGF 
calculations. On the other hand, to represent the electronic continuum within the timescale 
of simulation, a sufficient number of electronic states are required. As shown in Fig. [TlT c). 
20 states per lead is inadequate for both the transient I{t) or the steady state current. For 
the case of Fig. [TT] convergence is achieved when the number of states per lead is greater 
than 40. 



V. CONCLUDING REMARKS 



In this paper we have employed the ML-MCTDH-SQR method to simulate vibra- 
tionally coupled electron transport through single-molecule junctions. The ML-MCTDH- 
SQR method allows an accurate, in principle numerically exact treatment of this many-body 
quantum transport problem. The results obtained for a generic model demonstrate the im- 
portance of electronic- vibrational coupling, which has a significant infiuence on the transport 
properties. For situations where the energy of the bridge state is located close to the Fermi 
energy, the simulations show the time-dependent formation of a polaron state that results in 
a pronounced suppression of the current corresponding to the phenomenon of phonon block- 
ade. We have shown that this phenomenon cannot be explained solely by the polaron shift 
of the energy but requires methods that incorporate the dynamical effect of the vibrations 
on the transport. 

While some of these results have been discussed previously based on approximate meth- 
ods, the present methodology does not involve any systematic approximations and provides 
accurate benchmark results. It can thus also been used to test the validity of more approxi- 
mate methods. As an example, we have have discussed the validity of a NEGF method in the 
parameter regime of phonon blockade, demonstrating that accurate methods such as ML- 
MCTDH arc necessary to study transport in the strong coupling regime. A more detailed 
study of the validity of approximate methods as well as the application of the methodology 
to investigate signatures of vibronic effects in experimental transport spectra, such as, e.g., 
inelastic electron tunneling spectroscopy, will be the subject of future work. 

It is also emphasized that the time-dependent treatment employed in the ML-MCTDH- 
SQR method provides not only information on the steady state but also on the transient 
dynamics and can thus also been used to study the infiuence of time-dependent electrical 
fields, such as, e.g., ac gate fields or optical pulses on the transport process. 

In the present study we have focused on the effect of electronic-vibrational coupling 
on transport in molecular junctions. Another important mechanism is electron-electron 
interaction. The extension of the ML-MCTDH-SQR to include explicit electron-electron 
interaction is currently under way. This may open the perspective to a comprehensive 
many-body treatment of nonequilibrium charge transport at the nanoscale. 
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FIG. 1: Time-dependent current I{t) at different temperatures. Tlie parameters are: Oe = 0.2eV, 
Pe = leV, Ed — Ef = 0.5eV, and V = 0.2V. Results in panel (a) have been obtained without 
coupling to the vibrations (A = 0). In panel (b) couplings to an Ohmic bath of vibrational modes 
with parameters are A = 2000cm~^ and coc = 500cm~^ is included. The red line in panel (a) shows 
the current for a purely electronic model (A = 0) as obtained from Landauer theory 
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FIG. 2: Time-dependent current I{t) for different coupling strengths to the vibrational bath as 
specified by the reorganization energy A at temperature T = and bias voltage V = 0.2V. All 
other parameters are the same as in Fig. . 
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FIG. 3: Time-dependent current I{t) for different coupling strengths to the vibrational bath as 
specified by the reorganization energy A at temperature T = and bias voltage V = O.IV. Except 
for Ed — Ef = — 0.5eV, all other parameters are the same as in Fig. mb). 
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FIG. 4: Current-voltage characteristics for vibrational coupled electron transport. The parameters 
for the electronic lead states are Oe = 0.2eV and /3e = leV. The vibrational parameters are 
A = 2000cm- 1 andcjc = SOOcm'^ The energy of the bridge state is located at: (a) E^—Ej = 0.5eV, 
(b) E^ — Ef = — 0.5eV. The dashed lines depict the current for a purely electronic model (A = 0) 
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FIG. 5: Dependence of the time-dependent electric current (a) and the population of the bridge 
state (b) on the initial state. The electronic parameters are cce = 0.2eV, /3e = leV, E^—Ef = 0.5eV, 
and V = O.IV. The vibrational parameters are A = 2000cm~i and ojc = 500cm-^ 
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FIG. 6: Time-dependent current l(t) in the phonon-blockade regime for different coupling strengths 
to the vibrational bath as specified by the reorganization energy A. The electronic parameters are 
= 0.2eV, /3e = leV, Ed — Ef = 0, and V = O.IV. The characteristic frequency of the vibrational 
bath is uJc = 500cm~^. 




FIG. 7: Dependence of the current-voltage characteristics in the phonon-blockade regime on the 
electronic- vibrational coupling strength. The electronic parameters are the same as in Fig. [6l with 
= 500cm^^. 
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FIG. 8: Current-voltage characteristics for vibrationally coupled electron transport in the phonon- 
blockade regime. The electronic parameters are the same as in Fig. [H The vibrational parameters 
are A = 2000cm^^ and c^c = 500cm^^. Shown are results of a purely electronic model employing 
the bare (full line) and the polaron-shifted (dashed line) energy of the discrete state, respectively, 
as well as results of a full vibrationally-coupled many-body ML-MCTDH-SQR calculation (dashed- 
dotted line). 
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FIG. 9: Comparison of results obtain with NEGF theory (dashed hnes) and the ML-MCTDH-SQR 
method (diamonds) for vibronic couphng parameters A = 500cm^^ (a), A = lOOOcm^^ (b), and 
A = 2000cm^^ (c). Ah other parameters are the same as in Fig. [6l In addition, results for a purely 
electronic model are shown (full lines). 
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FIG. 10: Dependence of the time-dependent current I{t) in the phonon-blockade regime on the 
characteristic frequency of the vibrational bath. The electronic parameters are ag = 0.2eV, f3e = 
leV, Ed — Ef = 0, and V = O.IV. The reorganization energy is A = 2000cm~^. 
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FIG. 11: Time-dependent current I{t) for the parameter set of Oe = 0.2eV, /3e = leV, — Ef = 
0.5eV, and V = O.IV. Convergence is shown with respect to: (a) the number of SPFs for the 
electronic SPs; (b) the number of bath modes; (c) the number of electronic states for each lead. 
Other variational parameters are described in the text. 



